/**
 * @file   project.h
 * @author ubuntu <dadi@ubuntu>
 * @date   Tue Apr 20 14:00:16 2021
 * 
 * @brief  Homework project01.
 * 
 * 
 */

#ifndef __DADI__PROJECTH__
#define __DADI__PROJECTH__

#include<iostream>
#include<string>
#include<map>
#include<vector>
#include<utility>
#include<cmath>
#include<lapacke.h>

class TimeIntergrator
{
protected:
    int p;//accuracy
    double k;//time-step
public:
    std::vector<double> RHSFunc(std::vector<double> u);
    std::vector<double> DiffRHSFunc(std::vector<double> u);//second order Taylor expansion
    virtual void solver(std::vector<double> initialvalue, int period) = 0;
};

std::vector<double> TimeIntergrator::RHSFunc(std::vector<double> u)
{
    std::vector<double> u_temp(6);
    double mu = 1.0 / 81.45;
    u_temp[0] = u[3];
    u_temp[1] = u[4];
    u_temp[2] = u[5];
    double t2 = u[0] + mu - 1.0;
    double t3 = u[1] * u[1] + u[2] * u[2] + t2 * t2;
    double t4 = sqrt(t3 * t3 * t3);
    double t5 = u[1] * u[1] + u[2] * u[2] + (t2 + 1.0) * (t2 + 1.0);
    double t6 = sqrt(t5 * t5 * t5);
    u_temp[3] = 2.0 * u[4] + u[0] - mu * t2 / t4 - (1.0 - mu) * (t2 + 1.0) / t6;
    u_temp[4] = -2.0 * u[3] + u[1] - mu * u[1] / t4 - (1.0 - mu) * u[1] / t6;
    u_temp[5] = - mu * u[2] / t4 - (1.0 - mu) * u[2] / t6;
    u = u_temp;
    return u;
};

std::vector<double> TimeIntergrator::DiffRHSFunc(std::vector<double> u)
{
    std::vector<double> u_temp(6);
    double mu = 1.0 / 81.45;
    u_temp[0] = u[3];
    u_temp[1] = u[4];
    u_temp[2] = u[5];
    double t2 = u[0] + mu - 1.0;
    double t3 = u[1] * u[1] + u[2] * u[2] + t2 * t2;
    double t4 = sqrt(t3 * t3 * t3);
    double t5 = u[1] * u[1] + u[2] * u[2] + (t2 + 1.0) * (t2 + 1.0);
    double t6 = sqrt(t5 * t5 * t5);
    u_temp[3] = 2.0 * u[4] + u[0] - mu * t2 / t4 - (1.0 - mu) * (t2 + 1.0) / t6;
    u_temp[4] = -2.0 * u[3] + u[1] - mu * u[1] / t4 - (1.0 - mu) * u[1] / t6;
    u_temp[5] = - mu * u[2] / t4 - (1.0 - mu) * u[2] / t6;
    std::vector<std::vector<double>> u_m;
    u_m.resize(6);
    for(int i = 0; i < 6; i++)
	u_m[i].resize(6);
    for(int i = 0; i < 6; i++)
	for(int j = 0; j < 6; j++)
	    u_m[i][j] = 0;
    u_m[0][3] = 1.0;
    u_m[1][4] = 1.0;
    u_m[2][5] = 1.0;
    double t7 = sqrt(t3*t3*t3*t3*t3);
    double t8 = sqrt(t5*t5*t5*t5*t5);
    u_m[3][0] = (mu-1.0)/t6 - mu/t4 + (3.0*mu*t2*t2)/t7 - (3.0*(mu-1.0)*(t2+1.0)*(t2+1.0))/t8 + 1.0;
    u_m[3][1] = 3.0*mu*u[1]*t2/t7 - 3.0*(mu-1.0)*u[1]*(t2+1.0)/t8;
    u_m[3][2] = 3.0*mu*u[2]*t2/t7 - 3.0*(mu-1.0)*u[2]*(t2+1.0)/t8;
    u_m[3][4] = 2.0;
    u_m[4][0] = 3.0*mu*u[1]*t2/t7 - 3.0*(mu-1.0)*u[1]*(t2+1.0)/t8;
    u_m[4][1] = (mu-1.0)/t6 - mu/t4 - 3.0*u[1]*u[1]*(mu-1.0)/t8 + 3.0*mu*u[1]*u[1]/t7 + 1.0;
    u_m[4][2] = 3.0*mu*u[1]*u[2]/t7 - 3.0*(mu-1.0)*u[1]*u[2]/t8;
    u_m[4][3] = -2.0;
    u_m[5][0] = 3.0*mu*u[2]*t2/t7 - 3.0*(mu-1.0)*u[2]*(t2+1.0)/t8;
    u_m[5][1] = 3.0*mu*u[1]*u[2]/t7 - 3.0*(mu-1.0)*u[1]*u[2]/t8;
    u_m[5][2] = (mu-1.0)/t6 - mu/t4 - 3.0*u[2]*u[2]*(mu-1.0)/t8 + 3.0*mu*u[2]*u[2]/t7;
    for(int i = 0; i < 6; i++)
    {
	u[i] = 0;
	for(int j = 0; j < 6; j++)
	    u [i] += u_m[i][j]*u_temp[j];
    }
    return u;
};

class Linear_multistep_method : public TimeIntergrator
{
protected:
    int s;//step
    std::vector<double> alpha;
    std::vector<double> beta;
public:
    virtual void set_parameter() = 0;
};

class Adams_Bashforth : public Linear_multistep_method
{
public:
    Adams_Bashforth();
    Adams_Bashforth(int _p);
    void set_parameter() override;
    void solver(std::vector<double> initialvalue, int period) override;
};

Adams_Bashforth::Adams_Bashforth()
{
    s = -1;
    p = -1;
    k = 0.0;
    alpha = {0.0};
    beta = {0.0};
};

Adams_Bashforth::Adams_Bashforth(int _p)
{
    s = -1;
    p = _p;
    k = 0.0;
    alpha = {0.0};
    beta = {0.0};
};

void Adams_Bashforth::set_parameter()
{
    if(p == 1)
    {
	s = 1;
	k = 1e-5;
	alpha = {-1.0};
	beta = {1.0};
    }
    else if(p == 2)
    {
	s = 2;
	k = 1e-5;
	alpha = {0.0, -1.0};
	beta = {-1.0/2.0, 3.0/2.0};
    }
    else if(p == 3)
    {
	s = 3;
	k = 1e-5;
	alpha = {0.0, 0.0, -1.0};
	beta = {5.0/12.0, -16.0/12.0, 23.0/12.0};
    }
    else if(p == 4)
    {
	s = 4;
	k = 1e-5;
	alpha = {0.0, 0.0, 0.0, -1.0};
	beta = {-9.0/24.0, 37.0/24.0, -59.0/24.0, 55.0/24.0};
    }
    else;
};

void Adams_Bashforth::solver(std::vector<double> initialvalue, int period)
{
    set_parameter();
    std::vector<std::vector<double> > u;
    u.resize(s+1);
    for(int i = 0; i < s+1; i++)
    {
	u[i].resize(6);
    }
    std::vector<double> f;
    std::vector<double> df;
    u[0] = initialvalue;
    f = RHSFunc(u[0]);
    df = DiffRHSFunc(u[0]);
    for(int i = 1; i < s; i++)
	for(int j = 0; j < 6; j++)
	    u[i][j] = u[0][j] + i*k*f[j] + i*i*k*k*df[j]/2.0;
    std::cout << "u = [" << std::endl;
    double t = 0.0;
    while(t < period)
    {
	for(int i = 0; i < s; i++)
	{
	    f = RHSFunc(u[i]);
	    for(int j = 0; j < 6; j++)
		u[s][j] = 0;
	    for(int j = 0; j < 6; j++)
		u[s][j] += k*beta[i]*f[j] - alpha[i]*u[i][j];
	}
	for(int i = 0; i < 6; i++)
    	    std::cout << u[s][i] << " ";
    	std::cout << ";" << std::endl;
    	for(int i = 0; i < s; i++)
    	    u[i] = u[i+1];
    	t = t + k;
    }
    std::cout << "];" << std::endl;

    //solution error
    /* std::cout << "u = ["; */
    /* for(int i = 0; i < 6; i++) */
    /* 	{ */
    /* 	    std::cout << u[s-1][i] - initialvalue[i] << " "; */
    /* 	} */
    /* std::cout << "];" << std::endl; */
};

class Adams_Moulton : public Linear_multistep_method
{
private:
    std::vector<std::vector<double>> Jacobi_matrix;
public:
    Adams_Moulton();
    Adams_Moulton(int _p);
    void set_parameter() override;
    void solver(std::vector<double> initialvalue, int period) override;
    std::vector<std::vector<double>> Jacobi(std::vector<double> u);
    std::vector<double> Newton_method(std::vector<double> u, std::vector<double> g);
};

Adams_Moulton::Adams_Moulton()
{
    s = -1;
    p = -1;
    k = 0.0;
    alpha = {0.0};
    beta = {0.0};
};

Adams_Moulton::Adams_Moulton(int _p)
{
    s = -1;
    p = _p;
    k = 0.0;
    alpha = {0.0};
    beta = {0.0};
};

void Adams_Moulton::set_parameter()
{
    if(p == 2)
    {
	s = 1;
	k = 1e-5;
	alpha = {-1.0};
	beta = {1.0/2.0, 1.0/2.0};
    }
    else if(p == 3)
    {
	s = 2;
	k = 1e-5;
	alpha = {0.0, -1.0};
	beta = {-1.0/12.0, 8.0/12.0, 5.0/12.0};
    }
    else if(p == 4)
    {
	s = 3;
	k = 1e-5;
	alpha = {0.0, 0.0, -1.0};
	beta = {1.0/24.0, -5.0/24.0, 19.0/24.0, 9.0/24.0};
    }
    else if(p == 5)
    {
	s = 4;
	k = 1e-5;
	alpha = {0.0, 0.0, 0.0, -1.0};
	beta = {-19.0/720.0, 106.0/720.0, -264.0/720.0, 646.0/720.0, 251.0/720.0};
    }
    else;
};

void Adams_Moulton::solver(std::vector<double> initialvalue, int period)
{
    set_parameter();
    std::vector<std::vector<double> > u;
    u.resize(s+1);
    for(int i = 0; i < s+1; i++)
    {
	u[i].resize(6);
    }
    std::vector<double> f;
    std::vector<double> df;
    std::vector<double> g(6);
    u[0] = initialvalue;
    f = RHSFunc(u[0]);
    df = DiffRHSFunc(u[0]);
    for(int i = 1; i < s; i++)
	for(int j = 0; j < 6; j++)
	    u[i][j] = u[0][j] + i*k*f[j] + i*i*k*k*df[j]/2.0;
    std::cout << "u = [" << std::endl;
    double t = 0.0;
    while(t < period)
    {
	for(int i = 0; i < s; i++)
	{
	    f = RHSFunc(u[i]);
	    for(int j = 0; j < 6; j++)
		g[j] = 0;
	    for(int j = 0; j < 6; j++)
		g[j] += k*beta[i]*f[j] - alpha[i]*u[i][j];
	}
	u[s] = Newton_method(u[s-1], g);
	for(int i = 0; i < 6; i++)
	    std::cout << u[s][i] << " ";
	std::cout << ";" << std::endl;
	for(int i = 0; i < s; i++)
	    u[i] = u[i+1];
	t = t + k;
    }
    std::cout << "];" << std::endl;

    /* //solution error */
    /* std::cout << "u = ["; */
    /* for(int i = 0; i < 6; i++) */
    /* 	{ */
    /* 	    std::cout << u[s-1][i] - initialvalue[i] << " "; */
    /* 	} */
    /* std::cout << "];" << std::endl; */
};

std::vector<std::vector<double>> Adams_Moulton::Jacobi(std::vector<double> u)
{
    double mu = 1.0 / 81.45;
    double t2 = u[0] + mu - 1.0;
    double t3 = u[1] * u[1] + u[2] * u[2] + t2 * t2;
    double t4 = sqrt(t3 * t3 * t3);
    double t5 = u[1] * u[1] + u[2] * u[2] + (t2 + 1.0) * (t2 + 1.0);
    double t6 = sqrt(t5 * t5 * t5);
    Jacobi_matrix.resize(6);
    for(int i = 0; i < 6; i++)
	Jacobi_matrix[i].resize(6);
    for(int i = 0; i < 6; i++)
	for(int j = 0; j < 6; j++)
	    Jacobi_matrix[i][j] = 0;
    Jacobi_matrix[0][0] = 1.0;
    Jacobi_matrix[0][3] = -k*beta[s];
    Jacobi_matrix[1][1] = 1.0;
    Jacobi_matrix[1][4] = -k*beta[s];
    Jacobi_matrix[2][2] = 1.0;
    Jacobi_matrix[2][5] = -k*beta[s];
    double t7 = sqrt(t3*t3*t3*t3*t3);
    double t8 = sqrt(t5*t5*t5*t5*t5);
    Jacobi_matrix[3][0] = -k*beta[s]*((mu-1.0)/t6 - mu/t4 + (3.0*mu*t2*t2)/t7 - (3.0*(mu-1.0)*(t2+1.0)*(t2+1.0))/t8 + 1.0);
    Jacobi_matrix[3][1] = -k*beta[s]*(3.0*mu*u[1]*t2/t7 - 3.0*(mu-1.0)*u[1]*(t2+1.0)/t8);
    Jacobi_matrix[3][2] = -k*beta[s]*(3.0*mu*u[2]*t2/t7 - 3.0*(mu-1.0)*u[2]*(t2+1.0)/t8);
    Jacobi_matrix[3][3] = 1.0;
    Jacobi_matrix[3][4] = -k*beta[s]*2.0;
    Jacobi_matrix[4][0] = -k*beta[s]*(3.0*mu*u[1]*t2/t7 - 3.0*(mu-1.0)*u[1]*(t2+1.0)/t8);
    Jacobi_matrix[4][1] = -k*beta[s]*((mu-1.0)/t6 - mu/t4 - 3.0*u[1]*u[1]*(mu-1.0)/t8 + 3.0*mu*u[1]*u[1]/t7 + 1.0);
    Jacobi_matrix[4][2] = -k*beta[s]*(3.0*mu*u[1]*u[2]/t7 - 3.0*(mu-1.0)*u[1]*u[2]/t8);
    Jacobi_matrix[4][3] = k*beta[s]*2.0;
    Jacobi_matrix[4][4] = 1.0;
    Jacobi_matrix[5][0] = -k*beta[s]*(3.0*mu*u[2]*t2/t7 - 3.0*(mu-1.0)*u[2]*(t2+1.0)/t8);
    Jacobi_matrix[5][1] = -k*beta[s]*(3.0*mu*u[1]*u[2]/t7 - 3.0*(mu-1.0)*u[1]*u[2]/t8);
    Jacobi_matrix[5][2] = -k*beta[s]*((mu-1.0)/t6 - mu/t4 - 3.0*u[2]*u[2]*(mu-1.0)/t8 + 3.0*mu*u[2]*u[2]/t7);
    Jacobi_matrix[5][5] = 1.0;
    return Jacobi_matrix;
};

std::vector<double> Adams_Moulton::Newton_method(std::vector<double> u, std::vector<double> g)
{
    double norm = 0;
    std::vector<double> newton = {0, 0, 0, 0, 0, 0};
    double F[6] = {0, 0, 0, 0, 0, 0};
    for(int i = 0; i < 6; i++)
	F[i] = -(u[i] - k*beta[s]*RHSFunc(u)[i] - g[i]);
    for(int i = 0; i < 6; i++)
	if(fabs(F[i]) > norm)
	    norm = fabs(F[i]);
	else;
    while(norm > 1e-8)
    {
	double J[36];
	for(int i = 0; i < 6; i++)
	    for(int j = 0; j < 6; j++)
		J[j*6+i] = Jacobi(u)[i][j];
	int ipiv[6];
	LAPACKE_dgesv(LAPACK_COL_MAJOR, 6, 1, J, 6, ipiv, F, 6);
	for(int i = 0; i < 6; i++)
	    newton[i] = u[i] + F[i];
	u = newton;
	for(int i = 0; i < 6; i++)
	    F[i] = -(u[i] - k*beta[s]*RHSFunc(u)[i] - g[i]);
	norm = 0;
	for(int i = 0; i < 6; i++)
	    if(fabs(F[i]) > norm)
		norm = fabs(F[i]);
	    else;
    }
    return u;
};

class BDF : public Linear_multistep_method
{
private:
    std::vector<std::vector<double>> Jacobi_matrix;
public:
    BDF();
    BDF(int _p);
    void set_parameter() override;
    void solver(std::vector<double> initialvalue, int period) override;
    std::vector<std::vector<double>> Jacobi(std::vector<double> u);
    std::vector<double> Newton_method(std::vector<double> u, std::vector<double> g);
};

BDF::BDF()
{
    s = -1;
    p = -1;
    k = 0.0;
    alpha = {0.0};
    beta = {0.0};
};

BDF::BDF(int _p)
{
    s = -1;
    p = _p;
    k = 0.0;
    alpha = {0.0};
    beta = {0.0};
};

void BDF::set_parameter()
{
    if(p == 1)
    {
	s = 1;
	k = 1e-5;
	alpha = {-1.0};
	beta = {0.0, 1.0};
    }
    else if(p == 2)
    {
	s = 2;
	k = 1e-5;
	alpha = {1.0/3.0, -4.0/3.0};
	beta = {0.0, 0.0, 2.0/3.0};
    }
    else if(p == 3)
    {
	s = 3;
	k = 1e-5;
	alpha = {-2.0/11.0, 9.0/11.0, -18.0/11.0};
	beta = {0.0, 0.0, 0.0, 6.0/11.0};
    }
    else if(p == 4)
    {
	s = 4;
	k = 1e-5;
	alpha = {3.0/25.0, -16.0/25.0, 36.0/25.0, -48.0/25.0};
	beta = {0.0, 0.0, 0.0, 0.0, 12.0/25.0};
    }
    else;
};

void BDF::solver(std::vector<double> initialvalue, int period)
{
    set_parameter();
    std::vector<std::vector<double> > u;
    u.resize(s+1);
    for(int i = 0; i < s+1; i++)
    {
	u[i].resize(6);
    }
    std::vector<double> f;
    std::vector<double> df;
    std::vector<double> g(6);
    u[0] = initialvalue;
    f = RHSFunc(u[0]);
    df = DiffRHSFunc(u[0]);
    for(int i = 1; i < s; i++)
	for(int j = 0; j < 6; j++)
	    u[i][j] = u[0][j] + i*k*f[j] + i*i*k*k*df[j]/2.0;
    std::cout << "u = [" << std::endl;
    double t = 0.0;
    while(t < period)
    {
	for(int i = 0; i < s; i++)
	{
	    f = RHSFunc(u[i]);
	    for(int j = 0; j < 6; j++)
		g[j] = 0;
	    for(int j = 0; j < 6; j++)
		g[j] += k*beta[i]*f[j] - alpha[i]*u[i][j];
	}
	u[s] = Newton_method(u[s-1], g);
	for(int i = 0; i < 6; i++)
	    std::cout << u[s][i] << " ";
	std::cout << ";" << std::endl;
	for(int i = 0; i < s; i++)
	    u[i] = u[i+1];
	t = t + k;
    }
    std::cout << "];" << std::endl;

    /* //solution error */
    /* std::cout << "u = ["; */
    /* for(int i = 0; i < 6; i++) */
    /* 	{ */
    /* 	    std::cout << u[s-1][i] - initialvalue[i] << " "; */
    /* 	} */
    /* std::cout << "];" << std::endl; */
};

std::vector<std::vector<double>> BDF::Jacobi(std::vector<double> u)
{
    double mu = 1.0 / 81.45;
    double t2 = u[0] + mu - 1.0;
    double t3 = u[1] * u[1] + u[2] * u[2] + t2 * t2;
    double t4 = sqrt(t3 * t3 * t3);
    double t5 = u[1] * u[1] + u[2] * u[2] + (t2 + 1.0) * (t2 + 1.0);
    double t6 = sqrt(t5 * t5 * t5);
    Jacobi_matrix.resize(6);
    for(int i = 0; i < 6; i++)
	Jacobi_matrix[i].resize(6);
    for(int i = 0; i < 6; i++)
	for(int j = 0; j < 6; j++)
	    Jacobi_matrix[i][j] = 0;
    Jacobi_matrix[0][0] = 1.0;
    Jacobi_matrix[0][3] = -k*beta[s];
    Jacobi_matrix[1][1] = 1.0;
    Jacobi_matrix[1][4] = -k*beta[s];
    Jacobi_matrix[2][2] = 1.0;
    Jacobi_matrix[2][5] = -k*beta[s];
    double t7 = sqrt(t3*t3*t3*t3*t3);
    double t8 = sqrt(t5*t5*t5*t5*t5);
    Jacobi_matrix[3][0] = -k*beta[s]*((mu-1.0)/t6 - mu/t4 + (3.0*mu*t2*t2)/t7 - (3.0*(mu-1.0)*(t2+1.0)*(t2+1.0))/t8 + 1.0);
    Jacobi_matrix[3][1] = -k*beta[s]*(3.0*mu*u[1]*t2/t7 - 3.0*(mu-1.0)*u[1]*(t2+1.0)/t8);
    Jacobi_matrix[3][2] = -k*beta[s]*(3.0*mu*u[2]*t2/t7 - 3.0*(mu-1.0)*u[2]*(t2+1.0)/t8);
    Jacobi_matrix[3][3] = 1.0;
    Jacobi_matrix[3][4] = -k*beta[s]*2.0;
    Jacobi_matrix[4][0] = -k*beta[s]*(3.0*mu*u[1]*t2/t7 - 3.0*(mu-1.0)*u[1]*(t2+1.0)/t8);
    Jacobi_matrix[4][1] = -k*beta[s]*((mu-1.0)/t6 - mu/t4 - 3.0*u[1]*u[1]*(mu-1.0)/t8 + 3.0*mu*u[1]*u[1]/t7 + 1.0);
    Jacobi_matrix[4][2] = -k*beta[s]*(3.0*mu*u[1]*u[2]/t7 - 3.0*(mu-1.0)*u[1]*u[2]/t8);
    Jacobi_matrix[4][3] = k*beta[s]*2.0;
    Jacobi_matrix[4][4] = 1.0;
    Jacobi_matrix[5][0] = -k*beta[s]*(3.0*mu*u[2]*t2/t7 - 3.0*(mu-1.0)*u[2]*(t2+1.0)/t8);
    Jacobi_matrix[5][1] = -k*beta[s]*(3.0*mu*u[1]*u[2]/t7 - 3.0*(mu-1.0)*u[1]*u[2]/t8);
    Jacobi_matrix[5][2] = -k*beta[s]*((mu-1.0)/t6 - mu/t4 - 3.0*u[2]*u[2]*(mu-1.0)/t8 + 3.0*mu*u[2]*u[2]/t7);
    Jacobi_matrix[5][5] = 1.0;
    return Jacobi_matrix;
};

std::vector<double> BDF::Newton_method(std::vector<double> u, std::vector<double> g)
{
    double norm = 0;
    std::vector<double> newton = {0, 0, 0, 0, 0, 0};
    double F[6] = {0, 0, 0, 0, 0, 0};
    for(int i = 0; i < 6; i++)
	F[i] = -(u[i] - k*beta[s]*RHSFunc(u)[i] - g[i]);
    for(int i = 0; i < 6; i++)
	if(fabs(F[i]) > norm)
	    norm = fabs(F[i]);
	else;
    while(norm > 1e-8)
    {
	double J[36];
	for(int i = 0; i < 6; i++)
	    for(int j = 0; j < 6; j++)
		J[j*6+i] = Jacobi(u)[i][j];
	int ipiv[6];
	LAPACKE_dgesv(LAPACK_COL_MAJOR, 6, 1, J, 6, ipiv, F, 6);
	for(int i = 0; i < 6; i++)
	    newton[i] = u[i] + F[i];
	u = newton;
	for(int i = 0; i < 6; i++)
	    F[i] = -(u[i] - k*beta[s]*RHSFunc(u)[i] - g[i]);
	norm = 0;
	for(int i = 0; i < 6; i++)
	    if(fabs(F[i]) > norm)
		norm = fabs(F[i]);
	    else;
    }
    return u;
};

class Runge_Kutta_method : public TimeIntergrator
{
protected:
    std::vector<std::vector<double>> a;
    std::vector<double> b;
};

class classical_Runge_Kutta_method : public Runge_Kutta_method
{
public:
    classical_Runge_Kutta_method();
    void solver(std::vector<double> initialvalue, int period) override;
};

classical_Runge_Kutta_method::classical_Runge_Kutta_method()
{
    p = 4;
    k = 1e-5;
    a.resize(4);
    for(int i = 0; i < 4; i++)
	a[i].resize(4);
    for(int i = 0; i < 4; i++)
	for(int j = 0; j < 4; j++)
	    a[i][j] = 0;
    a[1][0] = 1.0/2.0;
    a[2][1] = 1.0/2.0;
    a[3][2] = 1.0;
    b = {1.0/6.0, 1.0/3.0, 1.0/3.0, 1.0/6.0};
};

void classical_Runge_Kutta_method::solver(std::vector<double> initialvalue, int period)
{
    std::vector<std::vector<double> > y;
    y.resize(4);
    for(int i = 0; i < 4; i++)
    {
	y[i].resize(6);
    }
    std::vector<double> u1(6);
    std::vector<double> u2(6);
    std::vector<double> temp(6);
    u1 = initialvalue;
    std::cout << "u = [" << std::endl;
    double t = 0.0;
    while(t < period)
    {
	y[0] = RHSFunc(u1);
	for(int i = 0; i < 6; i++)
	    temp[i] = u1[i] + k*y[0][i]*a[1][0];
	y[1] = RHSFunc(temp);
	for(int i = 0; i < 6; i++)
	    temp[i] = u1[i] + k*y[1][i]*a[2][1];
	y[2] = RHSFunc(temp);
	for(int i = 0; i < 6; i++)
	    temp[i] = u1[i] + k*y[2][i]*a[3][2];
	y[3] = RHSFunc(temp);
	for(int i = 0; i < 6; i++)
	    u2[i] = u1[i] + k*(b[0]*y[0][i]+b[1]*y[1][i]+b[2]*y[2][i]+b[3]*y[3][i]);
	for(int i = 0; i < 6; i++)
	    std::cout << u2[i] << " ";
	std::cout << ";" << std::endl;
	u1 = u2;
	t = t + k;
    }
    std::cout << "];" << std::endl;

    /* //solution error */
    /* std::cout << "u = ["; */
    /* for(int i = 0; i < 6; i++) */
    /* 	{ */
    /* 	    std::cout << u1[i] - initialvalue[i] << " "; */
    /* 	} */
    /* std::cout << "];" << std::endl; */
};


class TimeIntergratorFactory
{
public:
    typedef TimeIntergrator* (*TimeIntergratorCreator)();

public:
    bool RegisterTimeIntergrator(const std::pair<std::string, int> method_p, TimeIntergratorCreator createFn);
    bool UnregisterTimeIntergrator(const std::pair<std::string, int> method_p);
    TimeIntergrator* CreateTimeIntergrator(const std::pair<std::string, int> method_p) const;

private:
    typedef std::map<std::pair<std::string, int>, TimeIntergratorCreator> CallbackMap;
    CallbackMap callbacks_;

};

bool TimeIntergratorFactory::RegisterTimeIntergrator(const std::pair<std::string, int> method_p, TimeIntergratorCreator createFn)
{
    return callbacks_.insert(CallbackMap::value_type(method_p, createFn)).second;
};

bool TimeIntergratorFactory::UnregisterTimeIntergrator(const std::pair<std::string, int> method_p)
{
    return callbacks_.erase(method_p) == 1;
};

TimeIntergrator* TimeIntergratorFactory::CreateTimeIntergrator(const std::pair<std::string, int> method_p) const
{
// CallbackMap::const_iterator i = callbacks_.find(method_p);
    auto i = callbacks_.find(method_p);
    if (i == callbacks_.cend())
    { // not found
	throw std::runtime_error("Unknown TimeIntergrator ID");
    }
// Invoke the creation function
    return (i->second)();
};

TimeIntergrator* CreateAB1()
{
    return new Adams_Bashforth(1);
};
    
TimeIntergrator* CreateAB2()
{
    return new Adams_Bashforth(2);
};

TimeIntergrator* CreateAB3()
{
    return new Adams_Bashforth(3);
};

TimeIntergrator* CreateAB4()
{
    return new Adams_Bashforth(4);
};

TimeIntergrator* CreateAM2()
{
    return new Adams_Moulton(2);
};

TimeIntergrator* CreateAM3()
{
    return new Adams_Moulton(3);
};

TimeIntergrator* CreateAM4()
{
    return new Adams_Moulton(4);
};

TimeIntergrator* CreateAM5()
{
    return new Adams_Moulton(5);
};

TimeIntergrator* CreateBDF1()
{
    return new BDF(1);
};

TimeIntergrator* CreateBDF2()
{
    return new BDF(2);
};

TimeIntergrator* CreateBDF3()
{
    return new BDF(3);
};

TimeIntergrator* CreateBDF4()
{
    return new BDF(4);
};

TimeIntergrator* CreateRK4()
{
    return new classical_Runge_Kutta_method();
};

#else
//Do nothing
#endif
